Adaptive walks in a gene network model of morphogenesis: 
insights into the Cambrian explosion 



Ricard V. Sole/'^ Pau Fernandez/ and Stuart A. Kauffman^''' 

^ ICREA-Complex Systems Lab, Universitat Pompeu Fabra (GRIB), Dr Aiguader 80, 08003 Barcelona, Spain 
^Santa Fe Institute, 1399 Hyde Park Road, Santa Fe NM 87501, USA 

^Department of Cell Biology and Physiology, Health Science Center, University of New Mexico, 
Albuquerque, NM 87131 USA 



cn 

o 
o 

(N 
> 

o 

:z; 

o 



The emergence of complex patterns of organization close to the Cambrian boundary is known 
to have happened over a (geologically) short period of time. It involved the rapid diversification 
of body plans and stands as one of the major transitions in evolution. How it took place is a 
controversial issue. Here we explore this problem by considering a simple model of pattern forma- 
tion in multicellular organisms. By modeling gene network-based morphogenesis and its evolution 
through adaptive walks, we explore the question of how combinatorial explosions might have been 
actually involved in the Cambrian event. Here we show that a small amount of genetic complex- 
ity including both gene regulation and cell-cell signaling allows one to generate an extraordinary 
repertoire of stable spatial patterns of gene expression compatible with observed anteroposterior 
patterns in early development of metazoans. The consequences for the understanding of the tempo 
and mode of the Cambrian event are outlined. 



o 

d 
•1— I 

cr: 



> 



o 



I 



X 



Keywords: Eve devo; pattern formation; gene networks; morphogenesis; landscapes 



I. INTRODUCTION 

Morphological complexity in metazoa experienced an 
extraordinary leap close to the Cambrian boundary 
(around 550 Myrs ago). Such event has been labeled 
the Cambrian explosion. When this radiation began, 
and how rapi dly it u nfolde d, is still t he subject of ac- 
tive research (f Morris, "1998^ iRafil 11994: 'Sole et aILl2000t 
IValentine et aL. 1999) . The early origins of the last com- 
mon ancestor an d the structure of its bodyp lan are con- 
troversial issues l)Erwin and DavidsonLl200^ . The Cam- 
brian event established essentially all the major animal 
body plans and hence all the major phyla which would 
exist thereafter. A body plan can be described anatom- 
ically but also in terms of the spatiotempora l pattern 
of expression of some key genes llArthuilll99l . In this 
context, the origins of evolutionary novelty emerges as 
an integrative field inv olving gene regulation, dev elop- 
ment and paleobiology l|Arthuit l200l ICarrolll I200H) . As 
Shubin and Marshall pointed out, untangling the web 
of ecological, developmental and genetic interactions is a 
difficult task and a key question is which changes occur 
first ( Shubin and Marshall, 200Q.). 

A possible interpretation of the uniqueness, tempo and 
mode of the Cambrian event in terms of generic features 
of complex evolving systems was suggested by Kauffman 
in terms of adaptive wa lks on rugged fitness landscapes 
llKauffmanl.ll989lll99.1^ ■ In a nutshell, the idea is that, 
given the fundamental constraints imposed by early de- 
velopmental dynamics, early exploration of the universe 
of possible body plans took place quickly (once some un- 
derlying requirements were met) but then slowed down 
as the repertoire of possible bodyplans was filled out 
(Kauffman . 1989.; Raff , 19M)- Starting from some ini- 
tial condition where low-fit multicellular organisms were 
present, a rapid exploration of the landscape was allowed 



to occur. This initial exploration led to an increase of 
diversity of improved alternative morphologies thereby 
establishing phyla. As the rate of finding fitter mutants 
that alter early developmental processes (which define 
body plans) slowed down, lower taxonomic groups be- 
came established. 

The argument is completed by the assumption that 
mutants affecting early development have more profound 
effects than those affecting late development. In this sce- 
nario, by the Permian extinction (200 Myr later ), when 
an es timated 95% of species (82% of genera, see (jErwinl 
Il998|) ) went extinct, early developmental pathways would 
be expected to have become largely frozen in after the 
Cambrian event and no new phyla would be reachable. 
Fitter variants altering basic body plans would be very 
hard to find but not variants affecting late development. 
This would allow the radiation of new families and lower 
taxa. Fitness landscapes, first introduced by Wright in 
1932, allow to describe changes in the fitnes s of a given 
species through time in a well-defined fashion (Kau ffmanL 
E993.; Palme r, 1991; Rowe. .1994: .Stadler , 2002). In an 
idealized situation, we can imagine the landscape as a 
surface. Here the tops and bottoms of the landscape 
would indentify good and bad combinations of traits, re- 
spectively. This picture allows one to visualize the evo- 
lution of species as a hill climbing on the fitness surface. 

A more precise definition is provided by describing each 
individual (or species) is defined as a set of N binary 
variab les Sj € |0 , 1| defining a set of characteristic traits 
( K auffmanl. Il993t iKauffman and Levin. .1987.) . The total 
number of combinations is 2^ and each occupies a node 
of a iV-hypercube F^r. In figure Q we show an example 
of this fitness landscape for iV = 3. Each string has a 
well-defined fitness value, which can be represented by 
means of a TV-dimensional function (f> = (t>{Si, Sn)- If 
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FIG. 1 A simple fitness landscape witli A'^ = 3 traits and K — 2 epistatic interactions. Here the three traits involved interact 
with the other two (a). The possible states of the system are given by the vertices of the fitness cube (this time a A = 3- 
dimensional cube Fa). The resulting average fitness, <?!> = ^ (l>i/N is obtained from a fitness table (b) where the contribution of 
each trait <j)i (under the presence of the other two) is generated as a random number between zero and one. The table provides 
the local maxima. Adaptive walks (indicated as arrows) take place to nearest sites with higher fitness. In this example, adaptive 
walks end in two possible local peaks (c). 



the number of epistatic interactions K is zero, i. e. if dif- 
ferent traits do not influence each other, then a smooth 
landscape is obtained, whith a single peak. But if differ- 
ent traits interact, then the landscape becomes rugged 
and multiple local fitness peaks are allowed to occur. The 
simplest evolution on a fitness cube occurs by means sin- 
gle, one-bit steps. In other words, a given species can 
perform a random adaptive walk from a given node to- 
wards one of its N nearest neighbors if this leads to an 
increase in fitness (alternatively, a neutral change can 
also occur, and thus random drift is also allowed). A di- 
rect consequence of this assumption is that once a local 
peak is reached, no further changes are allowed. 

There is a universal featu re of adaptation on s t atisti- 
cally correlated landscapes (jKauffman and Levinl . Il987j) 
which can be appropriately used to test these ideas. This 
can be formulated as follows: if S is the cumulative num- 
ber of assumed improvements (new body plans) origi- 
nated and T is the cumulative number of tries, then S 
will grow with T in a logarithmic fashion. Graphically, 
this means that S grows rapidly at the beginning and 
then slows down. The cumulative number of tries can be 
approximated by the cumulative number of lower-level 
entities (e. g. genera) viewed as successive experiments 
in the generation of higher-level entities. The analysis of 
the cumulative increase of phylum originations against 
cumu lative genera seems to confirm this prediction (,,,Ebl3 . 
11998^ . 

A further step in exploring the possible scenarios for 
an explosion of patterns in the evolution of development 
should consider the developmental program in an explicit 
way. One important factor not addressed by previous 
theoretical studies deals with the spectrum of possible 
spatial patterns of gene expression that can be generated 
from a given number of genes. This is a relevant question 



since combinatorial explosions can occur once complex- 
ity thresholds (in number of genes or their interactions) 
are reached. In previous studies, it has been shown that 
some particular types of combinations of gene-gene inter- 
actions involving cell-cell communication allow to easily 
generate a number of spatial patterns including stripes , 
gradi ents or spots ijSalazar et all l200(l l200ll: ISole et gj] . 
I2OOI. These studies were performed on randomly wired 
networks and revealed a large fraction of spat ial patterns 
that could be generated and their nature l)Sole et all . 
|2002|) . One question immediately emerges: what are the 
minimal complexity requirements in terms of number of 
pattern-forming genes- in order to reach a high diversity 
of phenotypes? 

One possible approach to the problem is to con- 
sider an artificial evolution model were real organisms 
and their development ar e replaced by digital organisms 
l|Wilke and AdarnH l2002l) with simplified developmental 
programs. Such an approximation is becoming more 
common in evolutionary studies. The success of some of 
these models in extraordinary in providing insight into 
the evolution of complex biological systems. As pointed 
out by Gould, such a range of success is a consequence 
of universal l aws of change that are common to all com- 
plex systems (jGouldl f2003|) . In this context, we can con- 
jecture that digital developmental programs might also 
share generic properties with real, early development. 

Previous theoretical studies have shown that the fact 
of randomly wiring gene networks within a tissue context 
allows a high diversity of patterns and tha t spatial pat- 
terns are easily obtained ijSole et a^J . l2002j) . But no the- 
oretical study has addressed the question of the potential 
diversity of (stable) patterns that can be obtained by tun- 
ing the number of genes involved in creating them, pro- 
vided that some simple initial signal (an activated gene) 
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FIG. 2 (A) Gene interactions in early development. These 
is a subset of gene regulatory interactions that take place in 
Drosophila development. Arrows indicate diflterent types of 
positive and negative interactions among different genes. In 
(B) the simplified, discrete threshold model considered here 
is shown. Each gene is treated as a binary (Boolean) variable 
with only two states: active (1) or inactive (0). It integrates 
the input signals which are weighted through a matrix of links 
W (either positive or negative) . The output Y is obtained by 
using a threshold function, such as the one indicated in (C). 



is present at the beginning of development. Here we will 
explore this question by means of an simplified model of 
early development in which the number of different cell 
types is optimized. In this context, we specifically ask 
the following questions: 

1. What are the consequences of searching for stable 
patterns of increasing complexity, being complexity 
measured in terms of the number of cell types? 

2. What are the minimal requirements in terms of 
gene regulation complexity in order to be able to 
reach a wide spectrum of spatial patterns? 

3. What type of spatial patterns are obtained? 
4 



5. 



How does the evolution of these patterns takes 
place in terms of the underlying fitness landscape? 

Can a combinatorial explosion partially explain the 
tempo and mode of the Cambrian event? 



As will be shown below, the model approach provides 
nontrivial, tentative answers to the previous questions. 



II. GENE NETWORK MODEL 

A complete model of the gene activity pattern even 
at early stages of development would require a detailed 
descriptio n of the different levels of ge ne regulation and 
signaling (jArnone and Davidsonl Il997j) . Such a descrip- 
tion should take into account the continuous nature of 
mRNA and protein levels as well as other considerations 




FIG. 3 Modeling morphogenesis using simple gene network 
models. Here the simulated organism is a one-dimensional 
array of cells. Each cell contains the same set of A'^ genes, G 
of which (indicated by blue balls) only interact inside the cell 
whereas H others (red balls) act as microhormones, exchang- 
ing information with neighboring cells. 



related to the nature and distribution of cell signaling 
molecules and the stochastic nature of their dynamics. 
In this context, abstract models only taking into account 
a small amount of key features often captu re the essen- 
tial ingredients of gene regulation dyn amics ijHastv et all 
l200ll:ISole and Pastor-Satorrasll200^ . 

As a matter of fact, it has been recently shown 
that discrete, ON-OFF models of early development 
can actually provide complete enough information in 
order to re produce the key tra i ts of a developmen- 
tal pattern dAl bert and Othmcr, 2003: Bodnar. UOOl 
iMendoza et oITiigQa: .Sanchez and Thieffrv. .2001) . As 
discussed in l)Albert and Othmerl l2003|r within the con- 
text of the expression pattern of segment polarity genes 
in Drosophila, the observed patterns are determined by 
the topology of the network and the type of regulatory 
interactions between components. This is consistent with 
recent computer models indicating that a robust segment 
polarity module exists and that it is rather insensi tive to 
variation i n kinetic parameters llvon Dassow et 
(see also l|GibsonL l2002t iMeir et all |2002() ). In other 
words, a complete description of both wild type patterns 
and various mutants is successfully reached by using sim- 
ple Boolean networks. Here we take the same approach. 

Several choices can be made in dealing with the struc- 
ture of the wiring scheme to be used in the regulation 
network. One possibility is to restrict ourselves to a hi- 
erarchical cascade of interactions inspired in the topol- 
ogy of interactions of Drosophila. But it is kn own that 
even gap genes are not completely hierarchical ijGehrinel 
Il998|) . Some of them may act to switch on others (fig- 
ure 2a). Actually, this property characterizes other gene 
groups as well, reminding us that we are actually dealing 
with a complex network instead of a linear chain of steps 
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llArthud.lT997[) . 

The model explored here is a gene threshold system, 
described by a set of N genes per cell, interacting through 
a one-dimensional domain involving C cells, as shown in 
figure 13 The set of genes involved represent those re- 
quired to build the ba sic structure and are thus bauplan 
genes (as defined in l lTautz and Schmidl . I1998M Here 
gene states will be Boolean: genes are either active (1) 
or inactive (0) . Two types of elements will be considered: 
G genes and H hormones, with N = G + H. Genes inter- 
act within the cell, whose state at a given time t will be 
indicated as (t), with i — 1, G as the gene number 
and j = 1,...,C as the cell number (ordered from an- 
terior to posterior). Generically labeled microhormones 
l|,Tackson et oil Il98 6^ involve some implicit local media- 
tors communicating neighboring cells. These hormones 
can receive inputs from any of the first G units, but they 
can only make output to genes in other cells. The state of 
these H hormones will be accordingly indicated as (t) . 

Two matrices will be required in order to define 
the whole spectrum of links between different elements. 
These two matrices will be indicated by A = (Aik) and 
B = (i?ifc), defining interactions among the G genes and 
between genes and hormones, respectively. 

The basic set of equations of our gene network model 
are: 



giit + 1) = e[giit) + nm 
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gf(t) = Y^Akgfit) 
fc=i 



fc=i 



where S{x,y) is an "OR" function (that is, 5{x,y) = 1 
if either x = 1 or y = 1 and zero if x — y = 0). The 
function 8(a;) is a threshold function (that is, Q{x) — 1 
if X > and Q{x) — otherwise), and is shown in figure 

Therefore, genes are influenced either by genes or hor- 
mones and hormones are influenced just by genes, and 
the influences Qj and Tij add up to determine if genes 
will be active or inactive by the comparison to a threshold 
(in this case 0), as shown in figure 133. The OR function, 
or S{x, y), is used here as a substitute for the diffusion of 
microhormone to the neighbours, since a microhormonc 
will be seen as active in cell j if either of its neighbours 
is active. 

For cells at the boundaries we have the special terms: 



FIG. 4 The temporal construction of a pattern in a G = 3 
genes and H = 2 hormones network. (A) Scheme of represen- 
tation. (B) The 20 steps necessary to reach a stable pattern. 
Genes in black are active and genes in gray are inactive. 
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fe=l 



indicating that cells at the poles just have one neighbour. 

Finally, we need to define an initial condition. The 
simplest choice that can be defined involves the activation 
of a single gene at the anterior {j — 1) pole, as well as 
the symmetric case in which the signal is present at the 
two poles (j = 1 and j = C). Specifically, we set gl (0) = 
h{{0) = for all j = 1, ...,iVc, except ^^(0) = 1 in the 
single pole case, and an additional 5j^^(0) = 1 in the 
symmetric case. These choices correspond to maternal 
signals confined to the embryo's extremes. Such initial 
change will propagate to the rest of the tissue provided 
that the network is sensitive to it. 

As an example of the dynamics, the temporal evolution 
of one sample organism with one-pole maternal signal is 
shown in figure^ 
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III. ADAPTIVE WALKS 

Using the above mentioned gene network model of pat- 
tern formation, we now explore how pattern complexity 
emerges through a simple evolutionary algorithm where 
a population of P organisms evolve through adaptive 
walks. Such type of algorithm has been successfully 
used in different contexts (see for example (Niklas, 1994, 
[l£i3))- 

In order to define the evolution algorithm, we need to 
first define a fitness function. Here we restrict ourselves 
to measuring the number of different cell types. The 
number of cell types is a good measure of complexity 
which is known to increase through m etazoan evolution 
l|Carrol]l I2OOII: IValentine efaZI . Il994 . Increases in cell 
type number provide a high potential for further evolu- 
tion of anatomical and functional complexity, essentially 
through division of labor and the formation of spe cial- 
ized tissues ( Mavnard- Smith and Szathmarvl Il995|) . It 
is clear that the diversity of organisms involving a small 
number of cell types is fairly limited and thus a first step 
towards the evolution of complex organisms requires an 
expansion of their diversity of cell types. Actually, as dis- 
cussed by Erwin and Davidson, regulatory processes un- 
derlying cell-type specification are very old and display 
conserved plesiomorphic features ijErwin and DavidsonL 
|2flO,2). Morphogenetic processes would have evolved in- 
dependently and added at later stages. 

We consider, then, the number of cell types rict as our 
complexity measure and formulate the following ques- 
tion. What are the consequences for the diversity of spa- 
tial patterns that can be generated if we increase the 
number of cell types? 

Starting from a homogeneous population of P identi- 
cal organisms, we perform adaptive walks on the fitness 
landscape. At each generation in the algorithm, we se- 
quentially choose each organism and introduce a single 
change in its gene network. Three types of changes can 
occur, with probability ^, all of them affecting the de- 
pendencies between genes: 

• addition of new links, 

• removal of previously present links, and 

• randomization of links' weights. 

In this way the complexity of the network can be tuned 
independently for each organism. 

After the change, the developmental pattern generated 
by the organism is simulated. The number of different 
cell types n*^ is computed using this new pattern and it 
is compared with the previous one, If — "-c^^ 

the change is accepted and the new organism replaces 
the old one. Otherwise the change is rejected and the 
old organism is kept. The fact that we keep the new 
organism if n'j = 'T'**'^ (the exact developmental pat- 
tern might be different), introduces a certain amount of 
neutrality that may reduce the probability that an or- 
ganism gets trapped because it cannot find any changes 



that give it a higher number of cell types. As an ad- 
ditional requirement, the stability of the final pattern is 
enforced by rejecting those organisms which do not reach 
a stable state in / iterations. This is a strong constrain, 
since puts a hard limit on the attainable patterns, but it 
also seems sensible, for developmental processes are very 
robust and only display transient oscillations. 

The size of the space of possible patterns F is rather 
high, since the number of combinations of the activation 
of genes in each cell gives rise to a potential number of 
patterns of 

in = 2(^+«)^. (3) 

Actually, the tissue context exponentially expands the 
allowed space Tc of cell types reachable by isolated cells, 
which has a size \Tc\ = 2^"+'^) (i.e. we have |r| = \Taf). 
Even for a small number of cells and a small amount of 
genes, the resulting space is hyper astronomically large. 
For example, if we take G + H = A (which is one 
particularly interesting shown below) we have 

|r| = 2*^0 « 1017. 

In order to explore the potential repertoire of struc- 
tures that can be generated, we will focus in the set of 
stable spatial patterns Tg associated to each single gene 
(independently of others), from the set of possible spatial 
patterns, with jPgl = 2*^. This is equivalent to consider 
the raw capacity to generate specific patterns without 
regarding any gene as more important than any other. 
Therefore, a population of organisms at any given gen- 
eration has a corresponding set of patterns, V, the ones 
taken from all the patterns that the genes in each organ- 
ism generate individually in the developmental process. 
Any pattern in this list occurs in some specific genes in 
some organisms, across the different cells, and with a dif- 
ferent frequency, which we also measure. 

IV. RESULTS 

As an overall trend, simulations proceed in the direc- 
tion of increasing the number of patterns ["Pj, as the aver- 
age number of cell types grows generation by generation. 
Given that every organism in the population performs 
its own adaptive walk, one would expect the behavior of 
the average of the population would reflect some of the 
structure of the landscape. For instance, should there be 
any special local peaks in which an organism could get 
trapped, one would find many organisms there. As we 
will see, this is not the case. 

A graph of jT'l and the number of cell types achieved 
after 15000 generations for different numbers of genes and 
hormones is shown in figure [51 As the number of genes 
increases, the number of patterns discovered increases, 
with a tendency to saturate at high G. Since the total 
number of possible patterns is, for C = 9, |Fg| = 2^ = 
512, it is clear that the population almost finds all the 
possible patterns for high G. This is even more evident if 
we look at the number of cell types, which reaches a value 
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FIG. 5 (a) Number of stable spatial patterns obtained for 
different combinations of genes and microhormones from the 
adaptive walk algorithm. Here two sets of numerical experi- 
ments were performed, involving one and two hormones and 
varying numbers of genes. In (b) the corresponding numbers 
of cell types are shown for the two cases. Note the saturation 
towards the maximum number rict = C as the total number 
of genes increases. For each point, 5 simulations of 15000 gen- 
erations were performed with C = 9, P = 500, and / = 100. 



of almost 9, the maximum. From this curves it is also 
clear the simple change from one to two microhormones 
makes the system much less constrained and therefore 
more patterns are found. 

As mentioned earlier, given a population, many organ- 
isms generate some patterns more than once. In fact, the 
distribution F{r) of patterns generated has an interesting 
structure, as shown in figure Here patterns are sorted 
out by their rank r (i. e. by their order from the most 
to the least frequent). Some patterns, especially those 
with a small number of active cells at the anterior pole, 
are very frequent. Some others, on the contrary, are very 
scarce. The distribution follows approximately the form 
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FIG. 6 Frequency-rank distribution of complex patterns ob- 
tained from the adaptive walk evolutionary algorithm, for 
15000 generations, here with G = 3, = 2, P = 500 (a log-log 
plot is shown in the inset). A power decay is obtained, indi- 
cating a majority of patterns present in the final population 
but also a long tail of less common patterns. Some exam- 
ples are indicated. Most common patterns involve activation 
close to the anterior pole, but more complex patterns, such 
as stripes, appear to be rather frequent. As we move towards 
the tail of the distribution, more rich, asymmetric patterns 
are observed. 

of a power law, i. e. F{r) ^ r~", with a close to one. 

The frequent patters are, roughly, those that start the 
repeting process that generates a wave across the organ- 
ism but remain active only near the anterior pole. Since 
less combinations exist for that type of patterns, more or- 
ganisms will generate the same ones. The rare patterns 
seem to be those that are the result of the construction 
process and may vary considerably in its exact values, 
since many more different combinations exist. In figure 
^ the pattern generated by g2 would be an example of a 
"starter", and thus more frequent, and the gi would be 
an example of a possible rare one. 

Given a fixed G and H, two other parameters affect 
the number of patterns one obtains. On the one hand, 
the maximum number of iterations to attain a stable pat- 
tern, /, and, on the other, the number of organisms in the 
population, P. If / is smaller, less patterns are discov- 
ered since those patterns that may stabilize in a higher 
number of iterations are filtered out. If P is larger, the 
landscape is explored more extensively, and more pat- 
terns are discovered. Two figures demonstrate this de- 
pendence: figure[3shows two different runs with different 
population sizes all other parameters being equal, show- 
ing how the bigger population attains a bigger number 
of patterns; figure |H1 shows the exact dependency. 

The most surprising results, however, involve the 
repertoire of spatial patterns for small-sized gene net- 
works. The basic set of patterns obtained is shown in 
figure El for C = 15 cells, where the results from the two 
types of initial activations (one pole at left, two poles 
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FIG. 7 Time evolution of the number of stable patterns \V\ 
for two different population sizes. The other parameters are 
C = 15, G = 2, H — 2. The inset displays the same evolution 
in log-linear scale. A straight line in such a plot indicates 
a logarithmic increase in the number of patterns, consistent 
with a rugged landscape. 
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FIG. 8 The number of stable patterns \P\ for different popu- 
lation sizes. The other parameters are C = 15, G = 2, H = 2, 
and 15000 generations. The dashed line shows a curve that 
fits the data with a power function, i.e. f{x) — ax'^ with an 
exponent of /3 = 0.498 ± 0.02. 



at right) are shown. As we can see, the repertoire is 
fairly limited, even if we use G = 3 genes and H — 1 
hormones. Most patterns are rather homogeneous, dis- 
playing a small gradient or single stripes, although some 
patterns with wider stripes are also observable, thus indi- 
cating that the richness of dynamical patterns transiently 
generated by coupled genes can stably propagate even 
with only one hormone. 

The situation, however, dramatically changes once we 
cross a complexity threshold, li H = 2,G = 2, a very 
large number of patterns becomes available. In figure [TUI 
the resulting stable patterns are shown under the same 
simulation conditions. Now we can see that the num- 
ber of different patterns, |'Pi|=239 and 17^21 = 154 for one 



and two maternal signals, respectively, is comparable to 
the number of organisms involved. The diversity is also 
remarkable: many different orderings in the stripe and 
gap distributions are obtained, suggesting that a very 
high universe of combinations is compatible with a high 
diversity of cell types. 

V. DISCUSSION 

Any simple model of pattern formation does necessar- 
ily introduce shortcuts that limit the range of conclusions 
that can be safely reached. Our model does not include 
a large number of relevant features, from cell division to 
biologically realistic regulatory mechanisms. But again, 
as stressed in previous sections, some key, large scale 
features of complex systems do not depend on the spe- 
cific details involved at lower levels ijSole and Goodwinl 
120011) . This is fairly well exemplified by Boolean models 
of Drosophila development. 

Our results can be summarized as follows (providing a 
tentative list of answers to the list of questions presented 
in the first section): 

1. Using the number of cell types as a complexity 
measure, our model indicates that the whole spec- 
trum of spatial patterns is potentially reachable, 
provided that the population of digital organisms 
is large enough. Maximal diversity of cell types is 
positively correlated with the diversity of spatial 
patterns generated. This suggests that increasing 
cellular diversity is consistent with highly flexible 
pattern-forming mechanisms. 

2. If a single hormone is involved, complex patterns 
can be obtained by increasing the network com- 
plexity. A much more rapid increase is obtained by 
using two hormones. Actually, the combination of 
two genes -I- two hormones leads to a combinatorial 
explosion of spatial patterns. As a consequence, 
our results indicate that, provided the number of 
organisms is large enough, any pattern seems to be 
available. 

3. The spectrum of spatial patterns of gene expression 
is dominated by gap-like structures, stripes and all 
kind of combinations between them. It is remark- 
able that all classes of spatial structures are easily 
identified as matching the spatial distribution of 
maternal, gap and pair-rule gene expression pat- 
terns. Some genes (as it occurs in real develop- 
ment) are only transiently activated and thus ap- 
pear in the end as absent. 

4. The population climbs the underlying landscape in 
such a way that the rate of finding local peaks 
increases in a logarithmic fashion. This implies 
that a rapid diversification is followed by a further 
slowdown, in agreement with a rugged fitness land- 
scape. 
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FIG. 9 Taxonomy of patterns created from a system with different numbers of genes and microhormones at the subcritical 
regime. The left side shows the single pole case, whereas the right shows the symmetric one. A small amount of possible 
patterns is obtained, indicating that the possible repertoire of structures is fairly limited. Once two hormones are involved, it 
gets possible to obtain stable stripes, but only the most regular combinations are allowed. The experiment involved a population 
of P = 500 individuals for 15000 generations. 



5. The jump in pattern diversity experienced at H = 
G = 2 indicates that thresholds in network com- 
plexity, even at small-gene numbers exist and can 
lead to combinatorial explosions. Such explosions 
would open a whole spectrum of available struc- 
tures. Reaching such a threshold might require 
the formation of a minimal regulatory network and 
might also require other prerequisites dealing with 
body size, cellular interactions and tissue special- 
ization. However, once in place, the whole universe 
of patterns can be made suddenly available. 

The previous results support the idea that the Cam- 
brian event (and perhaps other rapid diversification 
events) might result from changes in the pattern of gene 
regulation. Of course the model does not consider ecolog- 
ical or other factors that might have played a key role. 



Instead, we concentrate in the study of the generative 
potential of such a simple model and in the universe of 
possible spatial patterns that can be generated. As a 
result, we obtain a surprising jump once a threshold of 
genetic complexity is reached. Although the whole reper- 
toire of patterns might need some exploratory effort to be 
reached, most patterns are easily found and once a large 
fraction of them has been obtained, further innovation 
occurs at a slower pace. 

The main message of this paper is that rapid diversifi- 
cation events in terms of generation of evolutionary nov- 
elty in developmental processes can take place through 
combinatorial explosions. Although it can be argued that 
the model is too simple, the exploration of continuous 
counterparts of these results give very similar outcomes 
(particularly the thesholds of network complexity). Ad- 
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FIG. 10 As in figure |^ but now using G = 2 and H = 2. A much higher diversity of patterns is generated, reveahng a 
combinatorial explosion of spatial motifs. The combination of complex patterns of gene interactions together with a combination 
of two microhormones leads to a great diversity of patterns of gene expression. The upper part shows the patterns for a single- 
pole maternal signal, whereas the lower part show the ones for the symmetric case. 
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ditionally, the continuous models allow to increase the 
repertoire of patterns and the same applies to discrete 
models using Boolean (instead of threshold) functions. 
In other words, the choices made here actually limit the 
diversity of patterns that can be generated. 
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